* SUBROUTINE MXDCMM               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX A
* BY A VECTOR X.
*
* PARAMETERS :
*  II  N  NUMBER OF ROWS OF THE MATRIX A.
*  II  M  NUMBER OF COLUMNS OF THE MATRIX A.
*  RI  A(N*M)  RECTANGULAR MATRIX STORED COLUMNWISE IN THE
*         ONE-DIMENSIONAL ARRAY.
*  RI  X(M)  INPUT VECTOR.
*  RO  Y(N)  OUTPUT VECTOR EQUAL TO A*X.
*
* SUBPROGRAMS USED :
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE MXDCMM(N,M,A,X,Y)
      INTEGER N,M
      DOUBLE PRECISION A(*),X(*),Y(*)
      INTEGER J,K
      CALL MXVSET(N,0.0D 0,Y)
      K=0
      DO 1 J=1,M
      CALL MXVDIR(N,X(J),A(K+1),Y,Y)
      K=K+N
    1 CONTINUE
      RETURN
      END
* SUBROUTINE MXDPGB                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A DENSE SYMMETRIC
* POSITIVE DEFINITE MATRIX A+E USING THE FACTORIZATION A+E=L*D*TRANS(L)
* OBTAINED BY THE SUBROUTINE MXDPGF.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
*         SUBROUTINE MXDPGF.
*  RU  X(N)  ON INPUT THE RIGHT HAND SIDE OF A SYSTEM OF LINEAR
*         EQUATIONS. ON OUTPUT THE SOLUTION OF A SYSTEM OF LINEAR
*         EQUATIONS.
*  II  JOB  OPTION. IF JOB=0 THEN X:=(A+E)**(-1)*X. IF JOB>0 THEN
*         X:=L**(-1)*X. IF JOB<0 THEN X:=TRANS(L)**(-1)*X.
*
* METHOD :
* BACK SUBSTITUTION
*
      SUBROUTINE MXDPGB(N,A,X,JOB)
      INTEGER JOB,N
      DOUBLE PRECISION A(*),X(*)
      INTEGER I,II,IJ,J
      IF (JOB.GE.0) THEN
*
*     PHASE 1 : X:=L**(-1)*X
*
          IJ = 0
          DO 20 I = 1,N
              DO 10 J = 1,I - 1
                  IJ = IJ + 1
                  X(I) = X(I) - A(IJ)*X(J)
   10         CONTINUE
              IJ = IJ + 1
   20     CONTINUE
      END IF
      IF (JOB.EQ.0) THEN
*
*     PHASE 2 : X:=D**(-1)*X
*
          II = 0
          DO 30 I = 1,N
              II = II + I
              X(I) = X(I)/A(II)
   30     CONTINUE
      END IF
      IF (JOB.LE.0) THEN
*
*     PHASE 3 : X:=TRANS(L)**(-1)*X
*
          II = N* (N-1)/2
          DO 50 I = N - 1,1,-1
              IJ = II
              DO 40 J = I + 1,N
                  IJ = IJ + J - 1
                  X(I) = X(I) - A(IJ)*X(J)
   40         CONTINUE
              II = II - I
   50     CONTINUE
      END IF
      RETURN
      END
* SUBROUTINE MXDPGD                ALL SYSTEMS                91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF A DIRECTION OF NEGATIVE CURVATURE WITH RESPECT TO A
* DENSE SYMMETRIC MATRIX A USING THE FACTORIZATION A+E=L*D*TRANS(L)
*         OBTAINED BY THE SUBROUTINE MXDPGF.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
*         SUBROUTINE MXDPGF.
*  RO  X(N)  COMPUTED DIRECTION OF NEGATIVE CURVATURE (I.E.
*         TRANS(X)*A*X<0) IF IT EXISTS.
*  II  INF  INFORMATION OBTAINED IN THE FACTORIZATION PROCESS. THE
*         DIRECTION OF NEGATIVE CURVATURE EXISTS ONLY IF INF>0.
*
* METHOD :
* P.E.GILL, W.MURRAY : NEWTON TYPE METHODS FOR UNCONSTRAINED AND
* LINEARLY CONSTRAINED OPTIMIZATION, MATH. PROGRAMMING 28 (1974)
* PP. 311-350.
*
      SUBROUTINE MXDPGD(N,A,X,INF)
      INTEGER N, INF
      DOUBLE PRECISION  A(N*(N+1)/2), X(N)
      INTEGER  I, J, II, IJ
      DOUBLE PRECISION ZERO, ONE
      PARAMETER (ZERO = 0.0D0, ONE = 1.0D0)
C
C     RIGHT HAND SIDE FORMATION
C
      DO 1 I=1,N
      X(I) = ZERO
    1 CONTINUE
      IF (INF .LE. 0) RETURN
      X(INF) = ONE
C
C     BACK SUBSTITUTION
C
      II = INF*(INF-1)/2
      DO  3  I = INF-1, 1, -1
      IJ = II
      DO  2  J = I+1, INF
      IJ = IJ + J - 1
      X(I) = X(I) - A(IJ)*X(J)
    2 CONTINUE
      II = II - I
    3 CONTINUE
      RETURN
      END
* SUBROUTINE MXDPGF                ALL SYSTEMS                89/12/01
* PORTABILITY : ALL SYSTEMS
* 89/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* FACTORIZATION A+E=L*D*TRANS(L) OF A DENSE SYMMETRIC POSITIVE DEFINITE
* MATRIX A+E WHERE D AND E ARE DIAGONAL POSITIVE DEFINITE MATRICES AND
* L IS A LOWER TRIANGULAR MATRIX. IF A IS SUFFICIENTLY POSITIVE
* DEFINITE THEN E=0.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RU  A(N*(N+1)/2)  ON INPUT A GIVEN DENSE SYMMETRIC (USUALLY POSITIVE
*         DEFINITE) MATRIX A STORED IN THE PACKED FORM. ON OUTPUT THE
*         COMPUTED FACTORIZATION A+E=L*D*TRANS(L).
*  IO  INF  AN INFORMATION OBTAINED IN THE FACTORIZATION PROCESS. IF
*         INF=0 THEN A IS SUFFICIENTLY POSITIVE DEFINITE AND E=0. IF
*         INF<0 THEN A IS NOT SUFFICIENTLY POSITIVE DEFINITE AND E>0. IF
*         INF>0 THEN A IS INDEFINITE AND INF IS AN INDEX OF THE
*         MOST NEGATIVE DIAGONAL ELEMENT USED IN THE FACTORIZATION
*         PROCESS.
*  RU  ALF  ON INPUT A DESIRED TOLERANCE FOR POSITIVE DEFINITENESS. ON
*         OUTPUT THE MOST NEGATIVE DIAGONAL ELEMENT USED IN THE
*         FACTORIZATION PROCESS (IF INF>0).
*  RO  TAU  MAXIMUM DIAGONAL ELEMENT OF THE MATRIX E.
*
* METHOD :
* P.E.GILL, W.MURRAY : NEWTON TYPE METHODS FOR UNCONSTRAINED AND
* LINEARLY CONSTRAINED OPTIMIZATION, MATH. PROGRAMMING 28 (1974)
* PP. 311-350.
*
      SUBROUTINE MXDPGF(N,A,INF,ALF,TAU)
      DOUBLE PRECISION ALF,TAU
      INTEGER INF,N
      DOUBLE PRECISION A(*)
      DOUBLE PRECISION BET,DEL,GAM,RHO,SIG,TOL
      INTEGER I,IJ,IK,J,K,KJ,KK,L
      L = 0
      INF = 0
      TOL = ALF
*
*     ESTIMATION OF THE MATRIX NORM
*
      ALF = 0.0D0
      BET = 0.0D0
      GAM = 0.0D0
      TAU = 0.0D0
      KK = 0
      DO 20 K = 1,N
          KK = KK + K
          BET = MAX(BET,ABS(A(KK)))
          KJ = KK
          DO 10 J = K + 1,N
              KJ = KJ + J - 1
              GAM = MAX(GAM,ABS(A(KJ)))
   10     CONTINUE
   20 CONTINUE
      BET = MAX(TOL,BET,GAM/N)
*      DEL = TOL*BET
      DEL = TOL*MAX(BET,1.0D0)
      KK = 0
      DO 60 K = 1,N
          KK = KK + K
*
*     DETERMINATION OF A DIAGONAL CORRECTION
*
          SIG = A(KK)
          IF (ALF.GT.SIG) THEN
              ALF = SIG
              L = K
          END IF
          GAM = 0.0D0
          KJ = KK
          DO 30 J = K + 1,N
              KJ = KJ + J - 1
              GAM = MAX(GAM,ABS(A(KJ)))
   30     CONTINUE
          GAM = GAM*GAM
          RHO = MAX(ABS(SIG),GAM/BET,DEL)
          IF (TAU.LT.RHO-SIG) THEN
              TAU = RHO - SIG
              INF = -1
          END IF
*
*     GAUSSIAN ELIMINATION
*
          A(KK) = RHO
          KJ = KK
          DO 50 J = K + 1,N
              KJ = KJ + J - 1
              GAM = A(KJ)
              A(KJ) = GAM/RHO
              IK = KK
              IJ = KJ
              DO 40 I = K + 1,J
                  IK = IK + I - 1
                  IJ = IJ + 1
                  A(IJ) = A(IJ) - A(IK)*GAM
   40         CONTINUE
   50     CONTINUE
   60 CONTINUE
      IF (L.GT.0 .AND. ABS(ALF).GT.DEL) INF = L
      RETURN
      END
* SUBROUTINE MXDPGN                ALL SYSTEMS               91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* ESTIMATION OF THE MINIMUM EIGENVALUE AND THE CORRESPONDING EIGENVECTOR
* OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX A+E USING THE
* FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE SUBROUTINE MXDPGF.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
*         SUBROUTINE MXDPGF.
*  RO  X(N)  ESTIMATED EIGENVECTOR.
*  RO  ALF  ESTIMATED EIGENVALUE.
*  II  JOB  OPTION. IF JOB=0 THEN ONLY ESTIMATED EIGENVALUE IS
*         COMPUTED. IS JOB>0 THEN BOTH ESTIMATED EIGENVALUE AND
*         ESTIMATED EIGENVECTOR ARE COMPUTED BY JOB ITERATIONS.
*
* SUBPROGRAMS USED :
*  S   MXDPGB  BACK SUBSTITUTION.
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*
* METHOD :
* A.K.CLINE, C.B.MOLER, G.W.STEWART, J.H.WILKINSON : AN ESTIMATE
* FOR THE CONDITION NUMBER OF A MATRIX. SIAM J. NUMER. ANAL. 16
* (1979) 368-373.
*
      SUBROUTINE MXDPGN(N,A,X,ALF,JOB)
      INTEGER N,JOB
      DOUBLE PRECISION A(N*(N+1)/2),X(N),ALF
      DOUBLE PRECISION XP,XM,FP,FM,MXVDOT
      INTEGER I,K,IK,KK
      DOUBLE PRECISION ZERO,ONE
      PARAMETER (ZERO=0.0D 0,ONE=1.0D 0)
C
C     COMPUTATION OF THE VECTOR V WITH POSSIBLE MAXIMUM NORM SUCH
C     THAT  L*D**(1/2)*V=U  WHERE U HAS ELEMENTS +1 OR -1
C
      DO 2 I=1,N
      X(I)=ZERO
    2 CONTINUE
      KK=0
      DO 6 K=1,N
      KK=KK+K
      XP=-X(K)+ONE
      XM=-X(K)-ONE
      FP=ABS(XP)
      FM=ABS(XM)
      IK=KK
      DO 3 I=K+1,N
      IK=IK+I-1
      FP=FP+ABS(X(I)+A(IK)*XP)
      FM=FM+ABS(X(I)+A(IK)*XM)
    3 CONTINUE
      IF(FP.GE.FM) THEN
      X(K)=XP
      IK=KK
      DO 4 I=K+1,N
      IK=IK+I-1
      X(I)=X(I)+A(IK)*XP
    4 CONTINUE
      ELSE
      X(K)=XM
      IK=KK
      DO 5 I=K+1,N
      IK=IK+I-1
      X(I)=X(I)+A(IK)*XM
    5 CONTINUE
      ENDIF
    6 CONTINUE
C
C     COMPUTATION OF THE VECTOR X SUCH THAT
C     D**(1/2)*TRANS(L)*X=V
C
      FM=ZERO
      KK=0
      DO 7 K=1,N
      KK=KK+K
      IF (JOB.LE.0) THEN
      FP=SQRT(A(KK))
      X(K)=X(K)/FP
      FM=FM+X(K)*X(K)
      ELSE
      X(K)=X(K)/A(KK)
      ENDIF
    7 CONTINUE
      FP=DBLE(N)
      IF (JOB.LE.0) THEN
C
C     FIRST ESTIMATION OF THE MINIMUM EIGENVALUE BY THE FORMULA
C     ALF=(TRANS(U)*U)/(TRANS(V)*V)
C
      ALF=FP/FM
      RETURN
      ENDIF
      FM=ZERO
      KK=N*(N+1)/2
      DO 9 K=N,1,-1
      IK=KK
      DO 8 I=K+1,N
      IK=IK+I-1
      X(K)=X(K)-A(IK)*X(I)
    8 CONTINUE
      FM=FM+X(K)*X(K)
      KK=KK-K
    9 CONTINUE
      FM=SQRT(FM)
      IF (JOB.LE.1) THEN
C
C     SECOND ESTIMATION OF THE MINIMUM EIGENVALUE BY THE FORMULA
C     ALF=SQRT(TRANS(U)*U)/SQRT(TRANS(X)*X)
C
      ALF=SQRT(FP)/FM
      ELSE
C
C     INVERSE ITERATIONS
C
      DO 11 K=2,JOB
C
C     SCALING THE VECTOR X BY ITS NORM
C
      DO 10 I=1,N
      X(I)=X(I)/FM
   10 CONTINUE
      CALL MXDPGB(N,A,X,0)
      FM=SQRT(MXVDOT(N,X,X))
   11 CONTINUE
      ALF=ONE/FM
      ENDIF
C
C     SCALING THE VECTOR X BY ITS NORM
C
      DO 12 I=1,N
      X(I)=X(I)/FM
   12 CONTINUE
      RETURN
      END
* FUNCTION MXDPGP                  ALL SYSTEMS                91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF THE NUMBER MXDPGP=TRANS(X)*D**(-1)*Y WHERE D IS A
* DIAGONAL MATRIX IN THE FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
* SUBROUTINE MXDPGF.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
*         SUBROUTINE MXDPGF.
*  RI  X(N)  INPUT VECTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RR  MXDPGP  COMPUTED NUMBER MXDPGP=TRANS(X)*D**(-1)*Y.
*
      FUNCTION MXDPGP(N,A,X,Y)
      INTEGER N
      DOUBLE PRECISION A(*),X(*),Y(*),MXDPGP
      DOUBLE PRECISION TEMP
      INTEGER I,J
      J = 0
      TEMP = 0.0D0
      DO 10 I = 1,N
          J = J + I
          TEMP = TEMP + X(I)*Y(I)/A(J)
   10 CONTINUE
      MXDPGP = TEMP
      RETURN
      END
* SUBROUTINE MXDPGS                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SCALING OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX A+E USING THE
* FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE SUBROUTINE MXDPGF.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RU  A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
*         SUBROUTINE MXDPGF.
*  RI  ALF  SCALING FACTOR.
*
      SUBROUTINE MXDPGS(N,A,ALF)
      DOUBLE PRECISION ALF
      INTEGER N
      DOUBLE PRECISION A(*)
      INTEGER I,J
      J = 0
      DO 10 I = 1,N
          J = J + I
          A(J) = A(J)*ALF
   10 CONTINUE
      RETURN
      END
* SUBROUTINE MXDPGU                ALL SYSTEMS                89/12/01
* PORTABILITY : ALL SYSTEMS
* 89/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* CORRECTION OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX A+E IN THE
* FACTORED FORM A+E=L*D*TRANS(L) OBTAINED BY THE SUBROUTINE MXDPGF.
* THE CORRECTION IS DEFINED AS A+E:=A+E+ALF*X*TRANS(X) WHERE ALF IS A
* GIVEN SCALING FACTOR AND X IS A GIVEN VECTOR.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RU  A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
*         SUBROUTINE MXDPGF.
*  RI  ALF  SCALING FACTOR IN THE CORRECTION TERM.
*  RI  X(N)  VECTOR IN THE CORRECTION TERM.
*  RA  Y(N) AUXILIARY VECTOR.
*
* METHOD :
* P.E.GILL, W.MURRAY, M.SAUNDERS: METHODS FOR COMPUTING AND MODIFYING
* THE LDV FACTORS OF A MATRIX, MATH. OF COMP. 29 (1974) PP. 1051-1077.
*
      SUBROUTINE MXDPGU(N,A,ALF,X,Y)
      DOUBLE PRECISION ZERO,ONE,FOUR,CON
      PARAMETER (ZERO=0.0D0,ONE=1.0D0,FOUR=4.0D0,CON=1.0D-8)
      DOUBLE PRECISION ALF,ALFR
      INTEGER N
      DOUBLE PRECISION A(*),X(*),Y(*)
      DOUBLE PRECISION B,D,P,R,T,TO
      INTEGER I,II,IJ,J
      IF (ALF.GE.ZERO) THEN
*
*     FORWARD CORRECTION IN CASE WHEN THE SCALING FACTOR IS NONNEGATIVE
*
          ALFR = SQRT(ALF)
          CALL MXVSCL(N,ALFR,X,Y)
          TO = ONE
          II = 0
          DO 30 I = 1,N
              II = II + I
              D = A(II)
              P = Y(I)
              T = TO + P*P/D
              R = TO/T
              A(II) = D/R
              B = P/ (D*T)
              IF (A(II).LE.FOUR*D) THEN
*
*     AN EASY FORMULA FOR LIMITED DIAGONAL ELEMENT
*
                  IJ = II
                  DO 10 J = I + 1,N
                      IJ = IJ + J - 1
                      D = A(IJ)
                      Y(J) = Y(J) - P*D
                      A(IJ) = D + B*Y(J)
   10             CONTINUE
              ELSE
*
*     A MORE COMPLICATE BUT NUMERICALLY STABLE FORMULA FOR UNLIMITED
*     DIAGONAL ELEMENT
*
                  IJ = II
                  DO 20 J = I + 1,N
                      IJ = IJ + J - 1
                      D = A(IJ)
                      A(IJ) = R*D + B*Y(J)
                      Y(J) = Y(J) - P*D
   20             CONTINUE
              END IF
              TO = T
   30     CONTINUE
      ELSE
*
*     BACKWARD CORRECTION IN CASE WHEN THE SCALING FACTOR IS NEGATIVE
*
          ALFR = SQRT(-ALF)
          CALL MXVSCL(N,ALFR,X,Y)
          TO = ONE
          IJ = 0
          DO 50 I = 1,N
              D = Y(I)
              DO 40 J = 1,I - 1
                  IJ = IJ + 1
                  D = D - A(IJ)*Y(J)
   40         CONTINUE
              Y(I) = D
              IJ = IJ + 1
              TO = TO - D*D/A(IJ)
   50     CONTINUE
          IF (TO.LE.ZERO) TO = CON
          II = N* (N+1)/2
          DO 70 I = N,1,-1
              D = A(II)
              P = Y(I)
              T = TO + P*P/D
              A(II) = D*TO/T
              B = -P/ (D*TO)
              TO = T
              IJ = II
              DO 60 J = I + 1,N
                  IJ = IJ + J - 1
                  D = A(IJ)
                  A(IJ) = D + B*Y(J)
                  Y(J) = Y(J) + P*D
   60         CONTINUE
              II = II - I
   70     CONTINUE
      END IF
      RETURN
      END
* SUBROUTINE MXDPRB                ALL SYSTEMS                89/12/01
* PORTABILITY : ALL SYSTEMS
* 89/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A DENSE SYMMETRIC
* POSITIVE DEFINITE MATRIX A USING THE FACTORIZATION A=TRANS(R)*R.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2) FACTORIZATION A=TRANS(R)*R.
*  RU  X(N)  ON INPUT THE RIGHT HAND SIDE OF A SYSTEM OF LINEAR
*         EQUATIONS. ON OUTPUT THE SOLUTION OF A SYSTEM OF LINEAR
*         EQUATIONS.
*  II  JOB  OPTION. IF JOB=0 THEN X:=A**(-1)*X. IF JOB>0 THEN
*         X:=TRANS(R)**(-1)*X. IF JOB<0 THEN X:=R**(-1)*X.
*
* METHOD :
* BACK SUBSTITUTION
*
      SUBROUTINE MXDPRB(N,A,X,JOB)
      INTEGER          JOB,N
      DOUBLE PRECISION A(*),X(*)
      INTEGER          I,II,IJ,J
      IF (JOB.GE.0) THEN
*
*     PHASE 1 : X:=TRANS(R)**(-1)*X
*
          IJ = 0
          DO 20 I = 1,N
              DO 10 J = 1,I - 1
                  IJ = IJ + 1
                  X(I) = X(I) - A(IJ)*X(J)
   10         CONTINUE
              IJ = IJ + 1
              X(I) = X(I)/A(IJ)
   20     CONTINUE
      END IF
      IF (JOB.LE.0) THEN
*
*     PHASE 2 : X:=R**(-1)*X
*
          II = N* (N+1)/2
          DO 40 I = N,1,-1
              IJ = II
              DO 30 J = I + 1,N
                  IJ = IJ + J - 1
                  X(I) = X(I) - A(IJ)*X(J)
   30         CONTINUE
              X(I) = X(I)/A(II)
              II = II - I
   40     CONTINUE
      END IF
      RETURN
      END
* SUBROUTINE MXDPRC                ALL SYSTEMS                92/12/01
* PORTABILITY : ALL SYSTEMS
* 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* CORRECTION OF A SINGULAR DENSE SYMMETRIC POSITIVE SEMIDEFINITE MATRIX
* A DECOMPOSED AS A=TRANS(R)*R.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RU  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM.
*  IO  INF  AN INFORMATION OBTAINED IN THE CORRECTION PROCESS. IF
*         INF=0 THEN A IS SUFFICIENTLY POSITIVE DEFINITE. IF
*         INF<0 THEN A IS NOT SUFFICIENTLY POSITIVE DEFINITE.
*         PROCESS.
*  RI  TOL  DESIRED TOLERANCE FOR POSITIVE DEFINITENESS.
*
      SUBROUTINE MXDPRC(N,A,INF,TOL)
      INTEGER N,INF
      DOUBLE PRECISION  A(N*(N+1)/2),TOL
      DOUBLE PRECISION TOL1,TEMP
      INTEGER L,I
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D 0)
      INF=0
      TOL1=SQRT(TOL)
      TEMP=TOL1
      DO 3 I=1,N*(N+1)/2
      TEMP=MAX(TEMP,ABS(A(I)))
    3 CONTINUE
      TEMP=TEMP*TOL1
      L=0
      DO 1 I=1,N
      L=L+I
      IF (ABS(A(L)).LE.TEMP) THEN
      A(L)=SIGN(TEMP,A(L))
      INF=-1
      ENDIF
    1 CONTINUE
      RETURN
      END
* SUBROUTINE MXDPRM                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF A GIVEN VECTOR X BY A DENSE SYMMETRIC POSITIVE
* DEFINITE MATRIX A USING THE FACTORIZATION A=TRANS(R)*R.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2) FACTORIZATION A=TRANS(R)*R.
*  RU  X(N)  ON INPUT THE GIVEN VECTOR. ON OUTPUT THE RESULT OF
*         MULTIPLICATION.
*  II  JOB  OPTION. IF JOB=0 THEN X:=A*X. IF JOB>0 THEN X:=R*X.
*         IF JOB<0 THEN X:=TRANS(R)*X.
*
      SUBROUTINE MXDPRM( N, A, X, JOB)
      INTEGER N, JOB
      DOUBLE PRECISION  A(N*(N+1)/2), X(N)
      INTEGER  I, J, II, IJ
      IF (JOB .GE. 0)  THEN
C
C     PHASE 1 : X:=R*X
C
      II = 0
      DO 3  I = 1, N
      II = II + I
      X(I) = A(II) * X(I)
      IJ = II
      DO 2  J = I+1, N
      IJ = IJ + J - 1
      X(I) = X(I) + A(IJ) * X(J)
    2 CONTINUE
    3 CONTINUE
      ENDIF
      IF (JOB .LE. 0)  THEN
C
C     PHASE 2 : X:=TRANS(R)*X
C
      IJ = N*(N+1)/2
      DO 6  I = N, 1, -1
      X(I) = A(IJ) * X(I)
      DO 5  J = I-1, 1, -1
      IJ = IJ - 1
      X(I) = X(I) + A(IJ) * X(J)
    5 CONTINUE
      IJ = IJ - 1
    6 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE MXDRGR               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* PLANE ROTATION IS APPLIED TO A ROWWISE STORED DENSE RECTANGULAR
* MATRIX A.
*
* PARAMETERS :
*  II  N  NUMBER OF COLUMNS OF THE MATRIX A.
*  RU  A(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE
*         ONE-DIMENSIONAL ARRAY.
*  II  K  FIRST INDEX OF THE PLANE ROTATION.
*  II  L  SECOND INDEX OF THE PLANE ROTATION.
*  RI  CK  DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX.
*  RI  CL  OFF-DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX.
*  II  IER  TYPE OF THE PLANE ROTATION. IER=0-GENERAL PLANE ROTATION.
*         IER=1-PERMUTATION. IER=2-TRANSFORMATION SUPPRESSED.
*
* SUBPROGRAMS USED :
*  S   MXVROT  PLANE ROTATION APPLIED TO TWO ELEMENTS.
*
      SUBROUTINE MXDRGR(N,A,K,L,CK,CL,IER)
      INTEGER N,K,L,IER
      DOUBLE PRECISION A(*),CK,CL
      INTEGER I,IK,IL
      IF (IER.NE.0.AND.IER.NE.1) RETURN
      IK=(K-1)*N
      IL=(L-1)*N
      DO 1 I=1,N
      IK=IK+1
      IL=IL+1
      CALL MXVROT(A(IK),A(IL),CK,CL,IER)
    1 CONTINUE
      RETURN
      END
* SUBROUTINE MXDRMD               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR MATRIX A BY
* A VECTOR X AND ADDITION OF A SCALED VECTOR ALF*Y.
*
* PARAMETERS :
*  II  N  NUMBER OF COLUMNS OF THE MATRIX A.
*  II  M  NUMBER OF ROWS OF THE MATRIX A.
*  RI  A(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE
*         ONE-DIMENSIONAL ARRAY.
*  RI  X(N)  INPUT VECTOR.
*  RI  ALF  SCALING FACTOR.
*  RI  Y(M)  INPUT VECTOR.
*  RO  Z(M)  OUTPUT VECTOR EQUAL TO A*X+ALF*Y.
*
      SUBROUTINE MXDRMD(N,M,A,X,ALF,Y,Z)
      INTEGER N,M
      DOUBLE PRECISION A(M*N),X(N),ALF,Y(M),Z(M)
      DOUBLE PRECISION TEMP
      INTEGER I,J,K
      K=0
      DO 2 J=1,M
      TEMP=ALF*Y(J)
      DO 1 I=1,N
      TEMP=TEMP+A(K+I)*X(I)
    1 CONTINUE
      Z(J)=TEMP
      K=K+N
    2 CONTINUE
      RETURN
      END
* SUBROUTINE MXDRMI               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* ROWWISE STORED DENSE RECTANGULAR MATRIX A IS SET TO BE A PART OF THE
* UNIT MATRIX.
*
* PARAMETERS :
*  II  N  NUMBER OF COLUMNS OF THE MATRIX A.
*  II  M  NUMBER OF ROWS OF THE MATRIX A.
*  RO  A(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE ONE-DIMENSIONAL
*          ARRAY. THIS MATRIX IS SET TO TRANS([I,0]).
*
      SUBROUTINE MXDRMI(N,M,A)
      INTEGER N,M
      DOUBLE PRECISION A(M*N)
      INTEGER I,J,K
      DOUBLE PRECISION ZERO,ONE
      PARAMETER (ZERO=0.0D 0,ONE=1.0D 0)
      K=0
      DO 2 J=1,M
      DO 1 I=1,N
      A(I+K)=ZERO
      IF (I.EQ.J) A(I+K)=ONE
    1 CONTINUE
      K=K+N
    2 CONTINUE
      RETURN
      END
* SUBROUTINE MXDRMM               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR MATRIX A BY
* A VECTOR X.
*
* PARAMETERS :
*  II  N  NUMBER OF COLUMNS OF THE MATRIX A.
*  II  M  NUMBER OF ROWS OF THE MATRIX A.
*  RI  A(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE
*         ONE-DIMENSIONAL ARRAY.
*  RI  X(N)  INPUT VECTOR.
*  RO  Y(M)  OUTPUT VECTOR EQUAL TO A*X.
*
      SUBROUTINE MXDRMM(N,M,A,X,Y)
      INTEGER N,M
      DOUBLE PRECISION A(*),X(*),Y(*)
      DOUBLE PRECISION TEMP
      INTEGER I,J,K
      K=0
      DO 2 J=1,M
      TEMP=0.0D 0
      DO 1 I=1,N
      TEMP=TEMP+A(K+I)*X(I)
    1 CONTINUE
      Y(J)=TEMP
      K=K+N
    2 CONTINUE
      RETURN
      END
* FUNCTION  MXDRMN               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* EUCLIDEAN NORM OF A PART OF THE I-TH COLUMN OF A ROWWISE STORED DENSE
* RECTANGULAR MATRIX A IS COMPUTED.
*
* PARAMETERS :
*  II  N  NUMBER OF COLUMNS OF THE MATRIX A.
*  II  M  NUMBER OF ROWS OF THE MATRIX A.
*  RI  A(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE
*         ONE-DIMENSIONAL ARRAY.
*  II  I  INDEX OF THE COLUMN WHOSE NORM IS COMPUTED.
*  II  J  INDEX OF THE FIRST ELEMENT FROM WHICH THE NORM IS COMPUTED.
*
      FUNCTION MXDRMN(N,M,A,I,J)
      INTEGER N,M,I,J
      DOUBLE PRECISION A(M*N),MXDRMN
      DOUBLE PRECISION POM,DEN
      INTEGER K,L
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D 0)
      DEN=ZERO
      L=(J-1)*N
      DO 1 K=J,M
      DEN=MAX(DEN,ABS(A(L+I)))
      L=L+N
    1 CONTINUE
      POM=ZERO
      IF (DEN.GT.ZERO) THEN
      L=(J-1)*N
      DO 2 K=J,M
      POM=POM+(A(L+I)/DEN)**2
      L=L+N
    2 CONTINUE
      ENDIF
      MXDRMN=DEN*SQRT(POM)
      RETURN
      END
* SUBROUTINE MXDRMV               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* K-TH COLUMN OF A ROWWISE STORED DENSE RECTANGULAR MATRIX A IS COPIED
* TO THE VECTOR X.
*
* PARAMETERS :
*  II  N  NUMBER OF COLUMNS OF THE MATRIX A.
*  II  M  NUMBER OF ROWS OF THE MATRIX A.
*  RI  A(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE
*         ONE-DIMENSIONAL ARRAY.
*  RO  X(M)  OUTPUT VECTOR SUCH THAT X(J)=A(J,K) FOR ALL J.
*  II  K  INDEX OF THE ROW BEING COPIED TO THE OUTPUT VECTOR.
*
      SUBROUTINE MXDRMV(N,M,A,X,K)
      INTEGER N,M,K
      DOUBLE PRECISION A(*),X(*)
      INTEGER I,J
      IF (K.LT.1.OR.K.GT.N) RETURN
      I=K
      DO 1 J=1,M
      X(J)=A(I)
      I=I+N
    1 CONTINUE
      RETURN
      END
* SUBROUTINE MXDRQF               ALL SYSTEMS                92/12/01
* PORTABILITY : ALL SYSTEMS
* 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* QR DECOMPOSITION OF ROWWISE STORED DENSE RECTANGULAR MATRIX Q USING
* HOUSEHOLDER TRANSFORMATIONS WITHOUT PIVOTING.
*
* PARAMETERS :
*  II  N  NUMBER OF COLUMNS OF THE MATRIX Q.
*  II  M  NUMBER OF ROWS OF THE MATRIX Q.
*  RU  Q(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE
*         ONE-DIMENSIONAL ARRAY.
*  RO  R(N*(N+1)/2)  UPPER TRIANGULAR MATRIX STORED IN THE PACKED FORM.
*
* SUBPROGRAMS USED :
*  S   MXDRMN  EUCLIDEAN NORM OF A PART OF THE ROWWISE STORED
*         RECTANGULAR MATRIX COLUMN.
*
* METHOD :
* P.A.BUSSINGER, G.H.GOLUB : LINEAR LEAST SQUARES SOLUTION BY
* HOUSEHOLDER TRANSFORMATION. NUMER. MATH. 7 (1965) 269-276.
*
      SUBROUTINE MXDRQF(N,M,Q,R)
      INTEGER N,M
      DOUBLE PRECISION Q(M*N),R(N*(N+1)/2)
      DOUBLE PRECISION ALF,POM,MXDRMN
      INTEGER I,J,K,L,JP,KP,NM
      DOUBLE PRECISION ZERO,ONE
      PARAMETER (ZERO=0.0D 0,ONE=1.0D 0)
      NM=MIN(N,M)
C
C     QR DECOMPOSITION
C
      L=0
      KP=0
      DO 6 K=1,NM
      POM=MXDRMN(N,M,Q,K,K)
      IF (Q(KP+K).NE.ZERO) POM=SIGN(POM,Q(KP+K))
      R(L+K)=-POM
      JP=0
      DO 1 J=1,K-1
      R(L+J)=Q(JP+K)
      Q(JP+K)=ZERO
      JP=JP+N
    1 CONTINUE
      IF (POM.NE.ZERO) THEN
C
C     HOUSEHOLDER TRANSFORMATION
C
      DO 2 J=K,M
      Q(JP+K)=Q(JP+K)/POM
      JP=JP+N
    2 CONTINUE
      Q(KP+K)=Q(KP+K)+ONE
      DO 5 I=K+1,N
      ALF=ZERO
      JP=KP
      DO 3 J=K,M
      ALF=ALF+Q(JP+K)*Q(JP+I)
      JP=JP+N
    3 CONTINUE
      ALF=ALF/Q(KP+K)
      JP=KP
      DO 4 J=K,M
      Q(JP+I)=Q(JP+I)-ALF*Q(JP+K)
      JP=JP+N
    4 CONTINUE
    5 CONTINUE
      ENDIF
      L=L+K
      KP=KP+N
    6 CONTINUE
C
C     EXPLICIT FORMULATION OF THE ORTHOGONAL MATRIX
C
      KP=N*N
      DO 11 K=N,1,-1
      KP=KP-N
      IF (Q(KP+K).NE.ZERO) THEN
      DO 9 I=K+1,N
      ALF=ZERO
      JP=KP
      DO 7 J=K,M
      ALF=ALF+Q(JP+K)*Q(JP+I)
      JP=JP+N
    7 CONTINUE
      ALF=ALF/Q(KP+K)
      JP=KP
      DO 8 J=K,M
      Q(JP+I)=Q(JP+I)-ALF*Q(JP+K)
      JP=JP+N
    8 CONTINUE
    9 CONTINUE
      JP=KP
      DO 10 J=K,M
      Q(JP+K)=-Q(JP+K)
      JP=JP+N
   10 CONTINUE
      ENDIF
      Q(KP+K)=Q(KP+K)+ONE
   11 CONTINUE
      RETURN
      END
* SUBROUTINE MXDRQU               ALL SYSTEMS                92/12/01
* PORTABILITY : ALL SYSTEMS
* 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* UPDATE OF A QR DECOMPOSITION. THIS QR DECOMPOSITION IS UPDATED
* BY THE RULE Q*R:=Q*R+ALF*X*TRANS(Y).
*
* PARAMETERS :
*  II  N  NUMBER OF COLUMNS OF THE MATRIX Q.
*  II  M  NUMBER OF ROWS OF THE MATRIX Q.
*  RU  Q(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE
*         ONE-DIMENSIONAL ARRAY (PART OF THE ORTHOGONAL MATRIX).
*  RU  R(N*(N+1)/2)  UPPER TRIANGULAR MATRIX STORED IN A PACKED FORM.
*  RI  ALF  SCALAR PARAMETER.
*  RI  X(M)  INPUT VECTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RA  Z(N)  AUXILIARY VECTOR.
*  IO  INF  INFORMATION. IF INF=0 THEN X LIES IN THE COLUMN SPACE OF Q.
*         IF INF=1 THEN X DOES NOT LIE IN THE COLUMN SPACE OF Q.
*
* SUBPROGRAMS USED :
*  RF  MXVNOR  EUCLIDEAN NORM OF A VECTOR.
*  S   MXVORT  COMPUTATION OF THE PLANE ROTATION MATRIX.
*  S   MXVROT  PLANE ROTATION IS APPLIED TO TWO NUMBERS.
*
* METHOD :
* J.W.DANIEL, W.B.GRAGG, L.KAUFMAN, G.W.STEWARD : REORTHOGONALIZATION
* AND STABLE ALGORITHMS FOR UPDATING THE GRAM-SCHMIDT QR FACTORIZATION.
* MATHEMATICS OF COMPUTATION 30 (1976) 772-795.
      SUBROUTINE MXDRQU(N,M,Q,R,ALF,X,Y,Z,INF)
      INTEGER N,M,INF
      DOUBLE PRECISION Q(M*N),R(N*(N+1)/2),ALF,X(M),Y(N),Z(N)
      DOUBLE PRECISION CK,CL,ZK,ZL,MXVNOR
      INTEGER J,K,L,KJ,KK,IER
      DOUBLE PRECISION ONE,CON
      PARAMETER (ONE=1.0D 0,CON=1.0D-10)
      INF=0
      KK=N*(N+1)/2
C
C     COMPUTATION OF THE VECTOR TRANS(Q)*X
C
      CALL MXDCMM(N,M,Q,X,Z)
      IF (M.GT.N) THEN
C
C     IF X DOES NOT LIE IN THE COLUMN SPACE OF Q WE HAVE TO USE
C     A SUBPROBLEM WHOSE DIMENSION IS BY ONE GREATER (INF=1).
C
      ZK=MXVNOR(M,X)
      CALL MXDRMD(N,M,Q,Z,-ONE,X,X)
      ZL=MXVNOR(M,X)
      IF (ZL.GT.CON*ZK) THEN
      INF=1
      CALL MXVSCL(M,-ONE/ZL,X,X)
      CALL MXVORT(Z(N),ZL,CK,CL,IER)
      IF (IER.EQ.0.OR.IER.EQ.1) THEN
      CALL MXVROT(R(KK),ZL,CK,CL,IER)
      KJ=N
      DO 1 J=1,M
      CALL MXVROT(Q(KJ),X(J),CK,CL,IER)
      KJ=KJ+N
    1 CONTINUE
      ENDIF
      ENDIF
      ENDIF
C
C     APPLICATION OF PLANE ROTATIONS TO THE VECTOR Z SO THAT
C     TRANS(Q1)*Z=E1 WHERE Q1 IS AN ORTHOGONAL MATRIX (ACCUMULATION OF
C     THE PLANE ROTATIONS) AND E1 IS THE FIRST COLUMN OF THE UNIT
C     MATRIX. AT THE SAME TIME BOTH THE UPPER HESSENBERG MATRIX
C     TRANS(Q1)*R AND THE ORTHOGONAL MATRIX Q*Q1 ARE CONSTRUCTED SO THAT
C     Q*Q1*R1=Q*Q1*(TRANS(Q1)*R+ALF*E1*TRANS(Y)) WHERE R1 IS AN UPPER
C     HESSENBERG MATRIX.
C
      DO 4 L=N,2,-1
      K=L-1
      KK=KK-L
      CALL MXVORT(Z(K),Z(L),CK,CL,IER)
      IF (IER.EQ.0.OR.IER.EQ.1) THEN
      CALL MXVROT(R(KK),Z(L),CK,CL,IER)
      KJ=KK
      DO 2 J=L,N
      KJ=KJ+J-1
      CALL MXVROT(R(KJ),R(KJ+1),CK,CL,IER)
    2 CONTINUE
      KJ=K
      DO 3 J=1,M
      CALL MXVROT(Q(KJ),Q(KJ+1),CK,CL,IER)
      KJ=KJ+N
    3 CONTINUE
      ENDIF
    4 CONTINUE
      Z(1)=ALF*Z(1)
      KJ=1
      DO 5 J=1,N
      R(KJ)=R(KJ)+Z(1)*Y(J)
      KJ=KJ+J
    5 CONTINUE
C
C     APPLICATION OF PLANE ROTATIONS TO THE UPPER HESSENBERG MATRIX R1
C     GIVEN ABOVE SO THAT R2=TRANS(Q2)*R1 WHERE Q2 IS AN ORTHOGONAL
C     MATRIX (ACCUMULATION OF THE PLANE ROTATIONS) AND R2 IS AN UPPER
C     TRIANGULAR MATRIX. WE OBTAIN THE NEW QR DECOMPOSITION Q*Q1*Q2*R2.
C
      KK=1
      DO 8 L=2,N
      K=L-1
      CALL MXVORT(R(KK),Z(L),CK,CL,IER)
      IF (IER.EQ.0.OR.IER.EQ.1) THEN
      KJ=KK
      DO 6 J=L,N
      KJ=KJ+J-1
      CALL MXVROT(R(KJ),R(KJ+1),CK,CL,IER)
    6 CONTINUE
      KJ=K
      DO 7 J=1,M
      CALL MXVROT(Q(KJ),Q(KJ+1),CK,CL,IER)
      KJ=KJ+N
    7 CONTINUE
      ENDIF
      KK=KK+L
    8 CONTINUE
C
C     BACK TRANSFORMATION OF THE GREATER SUBPROBLEM IF INF=1.
C
      IF (INF.EQ.1) THEN
      CALL MXVORT(R(KK),ZL,CK,CL,IER)
      IF (IER.EQ.0.OR.IER.EQ.1) THEN
      KJ=N
      DO 9 J=1,M
      CALL MXVROT(Q(KJ),X(J),CK,CL,IER)
      KJ=KJ+N
    9 CONTINUE
      ENDIF
      ENDIF
      RETURN
      END
* SUBROUTINE MXDSDA                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* A DENSE SYMMETRIC MATRIX A IS AUGMENTED BY THE SCALED UNIT MATRIX
* SUCH THAT A:=A+ALF*I (I IS THE UNIT MATRIX OF ORDER N).
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RU  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM.
*  RI  ALF  SCALING FACTOR.
*
      SUBROUTINE MXDSDA(N,A,ALF)
      INTEGER N
      DOUBLE PRECISION A(*),ALF
      INTEGER I,J
      J = 0
      DO  1  I = 1, N
      J = J + I
      A(J)=A(J)+ALF
    1 CONTINUE
      RETURN
      END
* FUNCTION MXDSDL                  ALL SYSTEMS                91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE MINIMUM DIAGONAL ELEMENT OF A DENSE SYMMETRIC
* MATRIX.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM.
*  IO  INF  INDEX OF THE MINIMUM DIAGONAL ELEMENT OF THE MATRIX A.
*  RR  MXDSDL  MINIMUM DIAGONAL ELEMENT OF THE MATRIX A.
*
      FUNCTION MXDSDL(N,A,INF)
      INTEGER N,INF
      DOUBLE PRECISION A(N*(N+1)/2),MXDSDL
      DOUBLE PRECISION TEMP
      INTEGER I,J
      J=1
      INF=1
      TEMP=A(1)
      DO 1 I=2,N
      J=J+I
      IF (TEMP.GT.A(J)) THEN
      INF=J
      TEMP=A(J)
      ENDIF
    1 CONTINUE
      MXDSDL=TEMP
      RETURN
      END
* SUBROUTINE MXDSMA             ALL SYSTEMS                 91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DENSE SYMMETRIC MATRIX AUGMENTED BY THE SCALED DENSE SYMMETRIC MATRIX.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRICES.
*  RI  ALF  SCALING FACTOR.
*  RI  A(N*(N+1)/2)  INPUT MATRIX.
*  RI  B(N*(N+1)/2)  INPUT MATRIX.
*  RO  C(N*(N+1)/2)  OUTPUT MATRIX WHERE C:=B+ALF*A.
*
      SUBROUTINE MXDSMA(N,ALF,A,B,C)
      INTEGER N
      DOUBLE PRECISION ALF,A(N*(N+1)/2),B(N*(N+1)/2),C(N*(N+1)/2)
      INTEGER I
      DO 1 I=1,N*(N+1)/2
      C(I)=B(I)+ALF*A(I)
    1 CONTINUE
      RETURN
      END
* SUBROUTINE MXDSMC                ALL SYSTEMS                91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COPYING OF A DENSE SYMMETRIC MATRIX.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRICES A AND B.
*  RI  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM.
*  RO  B(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM
*         WHERE B:=A.
*
      SUBROUTINE MXDSMC(N,A,B)
      INTEGER N
      DOUBLE PRECISION  A(N*(N+1)/2), B(N*(N+1)/2)
      INTEGER M,I
      M = N * (N+1)/2
      DO  1  I = 1, M
      B(I) = A(I)
    1 CONTINUE
      RETURN
      END
* SUBROUTINE MXDSMG                ALL SYSTEMS                91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
*  GERSHGORIN BOUNDS FOR EIGENVALUES OF A DENSE SYMMETRIC MATRIX.
*  AMIN .LE. ANY EIGENVALUE OF A .LE. AMAX.
*
* PARAMETERS :
*  II  N  DIMENSION OF THE MATRIX A.
*  RI  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM.
*  RO  AMIN  LOWER BOUND FOR EIGENVALUES OF A.
*  RO  AMAX  UPPER BOUND FOR EIGENVALUES OF A.
*
      SUBROUTINE MXDSMG(N,A,AMIN,AMAX)
      INTEGER  N
      DOUBLE PRECISION  A(N*(N+1)/2),AMIN,AMAX
      DOUBLE PRECISION TEMP
      INTEGER I,J,K,L
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D 0)
      AMAX=A(1)
      AMIN=A(1)
      K=0
      DO 3 I=1,N
      TEMP=ZERO
      L=K
      DO 1 J=1,I-1
      L=L+1
      TEMP=TEMP+ABS(A(L))
    1 CONTINUE
      L=L+1
      DO 2 J=I+1,N
      L=L+J-1
      TEMP=TEMP+ABS(A(L))
    2 CONTINUE
      K=K+I
      AMAX=MAX(AMAX,A(K)+TEMP)
      AMIN=MIN(AMIN,A(K)-TEMP)
    3 CONTINUE
      RETURN
      END
* SUBROUTINE MXDSMI                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DENSE SYMMETRIC MATRIX A IS SET TO THE UNIT MATRIX WITH THE SAME
* ORDER.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RO  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM
*         WHICH IS SET TO THE UNIT MATRIX (I.E. A:=I).
*
      SUBROUTINE MXDSMI(N,A)
      INTEGER          N
      DOUBLE PRECISION A(*)
      INTEGER          I,M
      M = N* (N+1)/2
      DO 10 I = 1,M
          A(I) = 0.0D0
   10 CONTINUE
      M = 0
      DO 20 I = 1,N
          M = M + I
          A(M) = 1.0D0
   20 CONTINUE
      RETURN
      END
* SUBROUTINE MXDSMM                ALL SYSTEMS                89/12/01
* PORTABILITY : ALL SYSTEMS
* 89/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF A DENSE SYMMETRIC MATRIX A BY A VECTOR X.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM.
*  RI  X(N)  INPUT VECTOR.
*  RO  Y(N)  OUTPUT VECTOR EQUAL TO  A*X.
*
      SUBROUTINE MXDSMM(N,A,X,Y)
      INTEGER          N
      DOUBLE PRECISION A(*),X(*),Y(*)
      DOUBLE PRECISION TEMP
      INTEGER          I,J,K,L
      K = 0
      DO 30 I = 1,N
          TEMP = 0.0D0
          L = K
          DO 10 J = 1,I
              L = L + 1
              TEMP = TEMP + A(L)*X(J)
   10     CONTINUE
          DO 20 J = I + 1,N
              L = L + J - 1
              TEMP = TEMP + A(L)*X(J)
   20     CONTINUE
          Y(I) = TEMP
          K = K + I
   30 CONTINUE
      RETURN
      END
* FUNCTION MXDSMQ                  ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VALUE OF A QUADRATIC FORM WITH A DENSE SYMMETRIC MATRIX A.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM.
*  RI  X(N)  GIVEN VECTOR.
*  RI  Y(N)  GIVEN VECTOR.
*  RR  MXDSMQ  VALUE OF THE QUADRATIC FORM MXDSMQ=TRANS(X)*A*Y.
*
      DOUBLE PRECISION
     +  FUNCTION MXDSMQ(N,A,X,Y)
      DOUBLE PRECISION ZERO
      PARAMETER        (ZERO=0.0D0)
      INTEGER          N
      DOUBLE PRECISION A(*),X(*),Y(*)
      DOUBLE PRECISION TEMP,TEMP1,TEMP2
      INTEGER          I,J,K
      TEMP = ZERO
      K = 0
      DO 20 I = 1,N
          TEMP1 = ZERO
          TEMP2 = ZERO
          DO 10 J = 1,I - 1
              K = K + 1
              TEMP1 = TEMP1 + A(K)*X(J)
              TEMP2 = TEMP2 + A(K)*Y(J)
   10     CONTINUE
          K = K + 1
          TEMP = TEMP + X(I)* (TEMP2+A(K)*Y(I)) + Y(I)*TEMP1
   20 CONTINUE
      MXDSMQ = TEMP
      RETURN
      END
* SUBROUTINE MXDSMR               ALL SYSTEMS                92/12/01
* PORTABILITY : ALL SYSTEMS
* 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* PLANE ROTATION IS APPLIED TO A DENSE SYMMETRIC MATRIX A. THE CASE
* K=L+1 IS REQUIRED.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RU  A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM.
*  II  K  FIRST INDEX OF PLANE ROTATION.
*  II  L  SECOND INDEX OF PLANE ROTATION.
*  RO  CK  DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX.
*  RO  CL  OFF-DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX.
*  IO  IER  INFORMATION ON THE TRANSFORMATION. IER<0-K OR L OUT OF
*         RANGE. IER=0-PLANE ROTATION. IER=1-PERMUTATION.
*         IER=2-TRANSFORMATION SUPPRESSED.
*
* SUBPROGRAMS USED :
*  S   MXVROT  PLANE ROTATION IS APPLIED TO TWO NUMBERS.
*
      SUBROUTINE MXDSMR(N,A,K,L,CK,CL,IER)
      INTEGER N,K,L,IER
      DOUBLE PRECISION A(*),CK,CL
      DOUBLE PRECISION AKK,AKL,ALL,CKK,CKL,CLL
      INTEGER J,KJ,LJ,KK,KL,LL
      IF (IER.NE.0.AND.IER.NE.1) RETURN
      IF (K.NE.L+1) THEN
      IER=-1
      RETURN
      ENDIF
      LJ=L*(L-1)/2
      DO 1 J=1,N
      IF (J.LE.L) THEN
      LJ=LJ+1
      KJ=LJ+L
      ELSE
      LJ=LJ+J-1
      KJ=LJ+1
      ENDIF
      IF (J.NE.K.AND.J.NE.L) THEN
      CALL MXVROT(A(KJ),A(LJ),CK,CL,IER)
      ENDIF
    1 CONTINUE
      IF (IER.EQ.0) THEN
      CKK=CK**2
      CKL=CK*CL
      CLL=CL**2
      LL=L*(L+1)/2
      KL=LL+L
      KK=LL+K
      AKL=(CKL+CKL)*A(KL)
      AKK=CKK*A(KK)+CLL*A(LL)+AKL
      ALL=CLL*A(KK)+CKK*A(LL)-AKL
      A(KL)=(CLL-CKK)*A(KL)+CKL*(A(KK)-A(LL))
      A(KK)=AKK
      A(LL)=ALL
      ELSE
      LL=L*(L+1)/2
      KK=LL+K
      AKK=A(KK)
      A(KK)=A(LL)
      A(LL)=AKK
      ENDIF
      RETURN
      END
* SUBROUTINE MXDSMS                ALL SYSTEMS                91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SCALING OF A DENSE SYMMETRIC MATRIX.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RU  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM
*         WHICH IS SCALED BY THE VALUE ALF (I.E. A:=ALF*A).
*  RI  ALF  SCALING FACTOR.
*
      SUBROUTINE MXDSMS( N, A, ALF)
      INTEGER N
      DOUBLE PRECISION  A(N*(N+1)/2), ALF
      INTEGER I,M
      M = N * (N+1) / 2
      DO 1  I = 1, M
      A(I) = A(I) * ALF
    1 CONTINUE
      RETURN
      END
* SUBROUTINE MXDSMU                ALL SYSTEMS                89/12/01
C PORTABILITY : ALL SYSTEMS
C 89/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* UPDATE OF A DENSE SYMMETRIC MATRIX A. THIS UPDATE IS DEFINED AS
* A:=A+ALF*X*TRANS(X) WHERE ALF IS A GIVEN SCALING FACTOR AND X IS
* A GIVEN VECTOR.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RU  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM.
*  RI  ALF  SCALING FACTOR IN THE CORRECTION TERM.
*  RI  X(N)  VECTOR IN THE CORRECTION TERM.
*
      SUBROUTINE MXDSMU( N, A, ALF, X)
      INTEGER N
      DOUBLE PRECISION  A(N*(N+1)/2), X(N), ALF
      DOUBLE PRECISION TEMP
      INTEGER  I, J, K
      K = 0
      DO 2  I = 1, N
      TEMP = ALF * X(I)
      DO 1  J = 1, I
      K = K + 1
      A(K) = A(K) + TEMP * X(J)
    1 CONTINUE
    2 CONTINUE
      RETURN
      END
* SUBROUTINE MXDSMV                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* K-TH ROW OF A DENSE SYMMETRIC MATRIX A IS COPIED TO THE VECTOR X.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM.
*  RO  X(N)  OUTPUT VECTOR.
*  II  K  INDEX OF COPIED ROW.
*
      SUBROUTINE MXDSMV(N,A,X,K)
      INTEGER          K,N
      DOUBLE PRECISION A(*),X(*)
      INTEGER          I,L
      L = K* (K-1)/2
      DO 10 I = 1,N
          IF (I.LE.K) THEN
              L = L + 1
          ELSE
              L = L + I - 1
          END IF
          X(I) = A(L)
   10 CONTINUE
      RETURN
      END
* SUBROUTINE MXVCOP                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COPYING OF A VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RO  Y(N)  OUTPUT VECTOR WHERE Y:= X.
*
      SUBROUTINE MXVCOP(N,X,Y)
      INTEGER N
      DOUBLE PRECISION X(*),Y(*)
      INTEGER I
      DO 10 I = 1,N
          Y(I) = X(I)
   10 CONTINUE
      RETURN
      END
* SUBROUTINE MXVDIF                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VECTOR DIFFERENCE.
*
* PARAMETERS :
*  RI  X(N)  INPUT VECTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RO  Z(N)  OUTPUT VECTOR WHERE Z:= X - Y.
*
      SUBROUTINE MXVDIF(N,X,Y,Z)
      INTEGER N
      DOUBLE PRECISION X(*),Y(*),Z(*)
      INTEGER I
      DO 10 I = 1,N
          Z(I) = X(I) - Y(I)
   10 CONTINUE
      RETURN
      END
* SUBROUTINE MXVDIR                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VECTOR AUGMENTED BY THE SCALED VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  A  SCALING FACTOR.
*  RI  X(N)  INPUT VECTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RO  Z(N)  OUTPUT VECTOR WHERE Z:= Y + A*X.
*
      SUBROUTINE MXVDIR(N,A,X,Y,Z)
      DOUBLE PRECISION A
      INTEGER N
      DOUBLE PRECISION X(*),Y(*),Z(*)
      INTEGER I
      DO 10 I = 1,N
          Z(I) = Y(I) + A*X(I)
   10 CONTINUE
      RETURN
      END
* FUNCTION MXVDOT                  ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DOT PRODUCT OF TWO VECTORS.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RR  MXVDOT  VALUE OF DOT PRODUCT MXVDOT=TRANS(X)*Y.
*
      FUNCTION MXVDOT(N,X,Y)
      INTEGER N
      DOUBLE PRECISION X(*),Y(*),MXVDOT
      DOUBLE PRECISION TEMP
      INTEGER I
      TEMP = 0.0D0
      DO 10 I = 1,N
          TEMP = TEMP + X(I)*Y(I)
   10 CONTINUE
      MXVDOT = TEMP
      RETURN
      END
* SUBROUTINE MXVINA             ALL SYSTEMS                   90/12/01
* PORTABILITY : ALL SYSTEMS
* 90/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* ELEMENTS OF THE INTEGER VECTOR ARE REPLACED BY THEIR ABSOLUTE VALUES.
*
* PARAMETERS :
*  II  N DIMENSION OF THE INTEGER VECTOR.
*  IU  IX(N)  INTEGER VECTOR WHICH IS UPDATED SO THAT IX(I):=ABS(IX(I))
*         FOR ALL I.
*
      SUBROUTINE MXVINA(N,IX)
      INTEGER N
      INTEGER IX(*)
      INTEGER I
      DO 10 I = 1,N
          IX(I) = ABS(IX(I))
          IF (IX(I).GT.10) IX(I) = IX(I) - 10
   10 CONTINUE
      RETURN
      END
* SUBROUTINE MXVIND               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* CHANGE OF THE INTEGER VECTOR ELEMENT FOR THE CONSTRAINT ADDITION.
*
* PARAMETERS :
*  IU  IX(N)  INTEGER VECTOR.
*  II  I  INDEX OF THE CHANGED ELEMENT.
*  II JOB  CHANGE SPECIFICATION. IS JOB.EQ.0 THEN IX(I)=10-IX(I).
*
      SUBROUTINE MXVIND(IX,I,JOB)
      INTEGER IX(*),I,JOB
      IF (JOB.EQ.0) IX(I)=10-IX(I)
      RETURN
      END
* SUBROUTINE MXVINS             ALL SYSTEMS                   90/12/01
C PORTABILITY : ALL SYSTEMS
C 90/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* INITIATION OF THE INTEGER VECTOR.
*
* PARAMETERS :
*  II  N DIMENSION OF THE INTEGER VECTOR.
*  II  IP  INTEGER PARAMETER.
*  IO  IX(N)  INTEGER VECTOR SUCH THAT IX(I)=IP FOR ALL I.
*
      SUBROUTINE MXVINS(N,IP,IX)
      INTEGER N,IP,IX(N)
      INTEGER I
      DO 1 I=1,N
      IX(I)=IP
    1 CONTINUE
      RETURN
      END
* SUBROUTINE MXVINV               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* CHANGE OF THE INTEGER VECTOR ELEMENT FOR THE CONSTRAINT ADDITION.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  IU  IX(N)  INTEGER VECTOR.
*  II  I  INDEX OF THE CHANGED ELEMENT.
*  II  JOB  CHANGE SPECIFICATION
*
      SUBROUTINE MXVINV(IX,I,JOB)
      INTEGER          I,JOB
      INTEGER          IX(*)
      IF ((IX(I).EQ.3.OR.IX(I).EQ.5) .AND. JOB.LT.0) IX(I) = IX(I) + 1
      IF ((IX(I).EQ.4.OR.IX(I).EQ.6) .AND. JOB.GT.0) IX(I) = IX(I) - 1
      IX(I) = -IX(I)
      RETURN
      END
* FUNCTION MXVMAX               ALL SYSTEMS                   91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* L-INFINITY NORM OF A VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RR  MXVMAX  L-INFINITY NORM OF THE VECTOR X.
*
      FUNCTION MXVMAX(N,X)
      INTEGER N
      DOUBLE PRECISION X(*),MXVMAX
      INTEGER I
      MXVMAX = 0.0D0
      DO 10 I = 1,N
          MXVMAX = MAX(MXVMAX,ABS(X(I)))
   10 CONTINUE
      RETURN
      END
* SUBROUTINE MXVNEG                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* CHANGE THE SIGNS OF VECTOR ELEMENTS.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RO  Y(N)  OUTPUT VECTOR WHERE Y:= - X.
*
      SUBROUTINE MXVNEG(N,X,Y)
      INTEGER N
      DOUBLE PRECISION X(*),Y(*)
      INTEGER I
      DO 10 I = 1,N
          Y(I) = -X(I)
   10 CONTINUE
      RETURN
      END
* FUNCTION  MXVNOR               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* EUCLIDEAN NORM OF A VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RR  MXVNOR  EUCLIDEAN NORM OF X.
*
      FUNCTION MXVNOR(N,X)
      INTEGER N
      DOUBLE PRECISION X(*),MXVNOR
      DOUBLE PRECISION POM,DEN
      INTEGER I
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D 0)
      DEN=ZERO
      DO 1 I=1,N
      DEN=MAX(DEN,ABS(X(I)))
    1 CONTINUE
      POM=ZERO
      IF (DEN.GT.ZERO) THEN
      DO 2 I=1,N
      POM=POM+(X(I)/DEN)**2
    2 CONTINUE
      ENDIF
      MXVNOR=DEN*SQRT(POM)
      RETURN
      END
* SUBROUTINE MXVORT               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF AN ELEMENTARY ORTHOGONAL MATRIX FOR PLANE ROTATION.
*
* PARAMETERS :
*  RU  XK  FIRST VALUE FOR PLANE ROTATION (XK IS TRANSFORMED TO
*         SQRT(XK**2+XL**2))
*  RU  XL  SECOND VALUE FOR PLANE ROTATION (XL IS TRANSFORMED TO
*         ZERO)
*  RO  CK  DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX.
*  RO  CL  OFF-DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX.
*  IO  IER  INFORMATION ON THE TRANSFORMATION. IER=0-GENERAL PLANE
*         ROTATION. IER=1-PERMUTATION. IER=2-TRANSFORMATION SUPPRESSED.
*
      SUBROUTINE MXVORT(XK,XL,CK,CL,IER)
      DOUBLE PRECISION CK,CL,XK,XL
      INTEGER          IER
      DOUBLE PRECISION DEN,POM
      IF (XL.EQ.0.0D0) THEN
          IER = 2
      ELSE IF (XK.EQ.0.0D0) THEN
          XK = XL
          XL = 0.0D0
          IER = 1
      ELSE
          IF (ABS(XK).GE.ABS(XL)) THEN
              POM = XL/XK
              DEN = SQRT(1.0D0+POM*POM)
              CK = 1.0D0/DEN
              CL = POM/DEN
              XK = XK*DEN
          ELSE
              POM = XK/XL
              DEN = SQRT(1.0D0+POM*POM)
              CL = 1.0D0/DEN
              CK = POM/DEN
              XK = XL*DEN
          END IF
          XL = 0.0D0
          IER = 0
      END IF
      RETURN
      END
* SUBROUTINE MXVROT               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* PLANE ROTATION IS APPLIED TO TWO VALUES.
*
* PARAMETERS :
*  RU  XK  FIRST VALUE FOR PLANE ROTATION.
*  RU  XL  SECOND VALUE FOR PLANE ROTATION.
*  RI  CK  DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX.
*  RI  CL  OFF-DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX.
*  II  IER  INFORMATION ON THE TRANSFORMATION. IER=0-GENERAL PLANE
*         ROTATION. IER=1-PERMUTATION. IER=2-TRANSFORMATION SUPPRESSED.
*
      SUBROUTINE MXVROT(XK,XL,CK,CL,IER)
      DOUBLE PRECISION CK,CL,XK,XL
      INTEGER          IER
      DOUBLE PRECISION YK,YL
      IF (IER.EQ.0) THEN
          YK = XK
          YL = XL
          XK = CK*YK + CL*YL
          XL = CL*YK - CK*YL
      ELSE IF (IER.EQ.1) THEN
          YK = XK
          XK = XL
          XL = YK
      END IF
      RETURN
      END
* SUBROUTINE MXVSAV                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DIFFERENCE OF TWO VECTORS RETURNED IN THE SUBSTRACTED ONE.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RU  Y(N)  UPDATE VECTOR WHERE Y:= X - Y.
*
      SUBROUTINE MXVSAV(N,X,Y)
      INTEGER          N
      DOUBLE PRECISION X(*),Y(*)
      DOUBLE PRECISION TEMP
      INTEGER          I
      DO 10 I = 1,N
          TEMP = Y(I)
          Y(I) = X(I) - Y(I)
          X(I) = TEMP
   10 CONTINUE
      RETURN
      END
* SUBROUTINE MXVSCL                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SCALING OF A VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RI  A  SCALING FACTOR.
*  RO  Y(N)  OUTPUT VECTOR WHERE Y:= A*X.
*
      SUBROUTINE MXVSCL(N,A,X,Y)
      DOUBLE PRECISION A
      INTEGER N
      DOUBLE PRECISION X(*),Y(*)
      INTEGER I
      DO 10 I = 1,N
          Y(I) = A*X(I)
   10 CONTINUE
      RETURN
      END
* SUBROUTINE MXVSET                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* A SCALAR IS SET TO ALL THE ELEMENTS OF A VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  A  INITIAL VALUE.
*  RO  X(N)  OUTPUT VECTOR SUCH THAT X(I)=A FOR ALL I.
*
      SUBROUTINE MXVSET(N,A,X)
      DOUBLE PRECISION A
      INTEGER N
      DOUBLE PRECISION X(*)
      INTEGER I
      DO 10 I = 1,N
          X(I) = A
   10 CONTINUE
      RETURN
      END
* SUBROUTINE MXVSUM                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SUM OF TWO VECTORS.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RO  Z(N)  OUTPUT VECTOR WHERE Z:= X + Y.
*
      SUBROUTINE MXVSUM(N,X,Y,Z)
      INTEGER          N
      DOUBLE PRECISION X(*),Y(*),Z(*)
      INTEGER          I
      DO 10 I = 1,N
          Z(I) = X(I) + Y(I)
   10 CONTINUE
      RETURN
      END
